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ABSTRACT 

Evolution of the mass function (MF) and radial distribution (RD) of the Galactic globular 
cluster (GC) system is calculated using an advanced and realistic Fokker-Planck (FP) model 
that considers dynamical friction, disc/bulge shocks, and eccentric cluster orbits. We perform 
hundreds of FP calculations with different initial cluster conditions, and then search a wide 
parameter space for the best-fitting initial GC MF and RD that evolves into the observed 
present-day Galactic GC MF and RD. By allowing both MF and RD of the initial GC system 
to vary, which is attempted for the first time in the present Letter, we find that our best-fitting 
models have a higher peak mass for a lognormal initial MF and a higher cut-off mass for a 
power-law initial MF than previous estimates, but our initial total masses in GCs, Mt,% = 1.5- 
1.8 x 10 8 Mq, are comparable to previous results. Significant findings include that our best- 
fitting lognormal MF shifts downward by 0.35 dex during the period of 13 Gyr, and that our 
power-law initial MF models well-fit the observed MF and RD only when the initial MF is 
truncated at > 10 5 M . We also find that our results are insensitive to the initial distribution 
of orbit eccentricity and inclination, but are rather sensitive to the initial concentration of the 
clusters and to how the initial tidal radius is defined. If the clusters are assumed to be formed 
at the apocentre while filling the tidal radius there, Mr,i can be as high as 6.9 x 10 s M , 
which amounts to ~ 75 per cent of the current mass in the stellar halo. 

Key words: stellar dynamics - Galaxy: globular clusters: general - Galaxy:evolution - 
Galaxy: formation - Galaxy: kinematics and dynamics 



1 INTRODUCTION 

Globular clusters (GCs) are the oldest (~ 13 Gyr) bound stellar 
subsystems in the Milky Way, and studies on the evolution of the 
GCs may give us valuable information on the environment of the 
Milky Way at the era of its formation. The initial mass function 
(MF) of the GC system is particularly of interest since it can tell 
us about the mode of star formation and about the fractions of stars 
that are formed inside and outside the clusters at the beginning of 
the galaxy. 

Evolution of the globular cluster mass function (GCMrO) is 
driven by many factors such as two-body relaxation, stellar evo- 
lution, binary heating, galactic tidal field, eccentric cluster orbits, 
and disc/bulge shocks. It is quite challenging to calculate the evo- 
lution of the GCMF considering all these factors. There have been 
numerous studies on the evolution of the GCMF using analytical 
models (Aguilar, Hut, & Ostriker 1988; Okazaki & Tosa 1995; 
Vesperini 1997; Fall & Zhang 2001, among others), Fokker-Planck 
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1 By the MF of the GC system or GCMF, we mean the number of clusters, 
not stars, as a function of mass. 



(FP) models (Gnedin & Ostriker 1997; Murali & Weinberg 1997, 
among others), and N-body models (Vesperini & Heggie 1997; 
Baumgardt 1998; Vesperini 1998, among others), but none of these 
studies took all of the aforementioned physics into account or im- 
plemented a wide enough range of parameter space for the initial 
conditions of the GCs. It is rather difficult to incorporate all these 
physics into analytical or FP models, whereas N-body simulations, 
although generally more accurate and easier to consider all the dis- 
ruption mechanisms than the former, are still too CPU-expensive to 
be performed for clusters with N > 10 , 

In the present letter, the evolution of the Galactic GCMFs is 
calculated using the most advanced, realistic FP model developed 
so far that incorporates all of the disruption mechanisms discussed 
earlier. We perform FP calculations for 720 different initial con- 
ditions (mass, galactocentric radius, orbit eccentricity, and orbit 
inclination). We then search a wide-parameter space for the best- 
fitting initial GC mass and radial distributions (RDs) that evolve 
into the observed present-day Galactic GC distributions (MF and 
RD). Such a simultaneous fit to both MF and RD of the GC system 
is attempted for the first time in the present study. We adopt a log- 
normal and a truncated power-law function for the initial GCMF 
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(IGCMF0 Our IGCMFs are to be regarded as the models after the 
gas expulsion., and a softened power-law function for the initial RD 
of the apocentre. 



2 THE PRESENT-DAY GCMF 

The Galactic GCs can be classified into three groups by their age 
and metallicity: the 'old' halo (OH) and bulge/disc (BD) clusters 
are believed to be Galactic natives that were created when the 
protogalaxy collapsed, while the 'young' halo (YH) clusters are 
thought to have been formed in external satellite galaxies (Zinn 
1993; Parmentier et al. 2000). When comparing our calculations 
to the present-day GCMF, we only consider the 'native' clusters, 
i.e., the OH and BD clusters. We adopt the GC classification by 
Mackey & van den Bergh (2005), which are based on the data base 
compiled by Harris (1996). Our native GCs do not include the six 
objects that belong to the Sagittarius dwarf, seven objects whose 
origins remain unknown (Mackey & van den Bergh 2005), and 15 
objects that are thought to be the remnants of dwarf galaxies (Lee, 
Gim, & Casetti-Dinescu 2007). The total number of the native clus- 
ters is 95 (61 OHs and 34 BDs), and the total mass in these clusters 
is 2.6 x 10 7 M Q (1.8 x 10 7 M Q in OHs and 0.8 x 10 7 M in BDs). 
When fit to a lognormal function, the MF of our native clusters has 
M p = 10 5 - 26 M Q and ciogM = 0.44. The observed luminosities 
of the GCs are transformed into masses with a mass-to-light ratio 
of M/L v = 2. 



3 MODELS AND INITIAL CONDITIONS 

We adopt the anisotropic FP model developed by Takahashi & Lee 
(2000) and Takahashi & Portegies Zwart (2000), which directly 
integrates the orbit- averaged FP equation of two (energy-angular 
momentum) dimensions and considers multiple stellar mass com- 
ponents and the effects of tidal fields, three-body binary heating and 
stellar evolution. To this model, we have added the effects of tidal 
binary heating, disc/bulge shocks, dynamical friction and realistic 
cluster orbits. 

Disk and bulge shocks arise when clusters pass through the 
galactic disc or the bulge. Shocks inject kinetic energy into the 
cluster and speedup its disruption. Gnedin, Lee, & Ostriker (1999) 
incorporated these shocks into an FP model of one (energy) di- 
mension. We have extended their recipe and applied it to our two- 
dimensional (2D) FP model (detailed description on this applica- 
tion will be presented in Shin, Kim, & Takahashi 2008). The FP 
model is numerically stable in most cases, but we find that it en- 
counters numerical problems rather often when the effects of tidal 
shocks are included in the anisotropic FP model. To avoid such a 
problem, Shin & Kim (2007) developed a new integration scheme 
for a 2D FP equation by adopting an Alternating Direction Implicit 
method. We use this scheme for our calculations. 

Dynamical friction between a cluster and galactic field stars 
gradually transports a cluster to the inner region of the galaxy. Since 
a cluster with a given mass M has a smaller tidal radius when lo- 
cated at a smaller galactocentric radius R, dynamical friction in- 
creases the mass loss rate M of the cluster over its tidal radius 

2 Parmentier & Gilmore (2007) find that a power-law initial mass distri- 
bution of protoglobular clouds can quickly evolve into a lognormal initial 
GCMF due to the expulsion of the gas remnant from star formation, if the 
power-law mass distribution has a lower mass limit 



r t . Eccentric cluster orbits have similar effects on the clusters, al- 
though the effects are transient and periodic. To incorporate the ef- 
fects of dynamical friction and eccentric orbits into our FP model, 
we follow the orbit of the cluster by integrating the equation of mo- 
tion in the Galactic potential with a drag due to dynamical friction, 
and continuously update r t of the cluster at each time-step, which is 
determined by the current 7? and M.0For a drag due to dynamical 
friction, we adopt the formula by Chandrasekhar (1943), and for 
the Galactic potential, we employ the model by Johnston, Spergel, 
&Hernquist(1995). 

Takahashi & Portegies Zwart (1998, 2000) found that FP mod- 
els produce results similar to those from N-body simulations at 
least for clusters on circular orbits, if an "apocenter criterion" and 
u eac = 2-3 are used for the escape criterion of the FP model 
(v esc is a dimensionless parameter that determines the time-scale 
on which escaping stars leave the cluster). We adopt the apocenter 
criterion as well and v asc = 2.5. We find that the FP and N-body 
methods show a good agreement for clusters on eccentric orbits as 
well. A comparison between our FP calculations and N-body sim- 
ulations of Baumgardt & Makino (2003) for clusters on eccentric 
orbits with initial masses larger than 10 4 Mq shows a good agree- 
ment of cluster lifetimes within ~ 25 per cent. 

Parameters for our FP survey are the following four initial 
cluster conditions: M, R, the inclination of the orbital plane i rel- 
ative to the Galactic plane, and the orbit eccentricity e, which is 
defined as (R a — R p )/(R a + R P ) with R a and R p being the apoc- 
entre and pericentre distances, respectively. An initial R is defined 
to be R a of the orbit (i.e., clusters are initially located at R a ). 

We choose eight M values from 10 s ' 5 to 10 7 Mq and six 
R values from 10 3 to 10 4 ' 67 pc, both equally spaced on the log- 
arithmic scale. For the orbit inclination and eccentricity, we choose 
e = 0, 0.125, 0.25, 0.5, and 0.75, and i = 15°, 45°, and 75°, 
respectively. We perform FP calculations for all possible combina- 
tions out of these four parameters, thus the total number of cluster 
models considered in our study amounts to 720. 

The stellar density and velocity dispersion distributions within 
each cluster follow the King model (King 1966) with a concentra- 
tion parameter Wo — 7 and with no initial velocity anisotropy and 
no mass segregation. Clusters are assumed to initially fill the tidal 
radius, but there is a question regarding what tidal radius one needs 
to adopt for clusters with eccentric orbits — if the star formation in 
a cluster takes place on a time-scale much longer than the orbital 
timescale of the cluster, the tidal radius at the pericentre would be 
appropriate, and if much shorter, any value between the pericen- 
tre and the apocentre would be possible. Here, we assume that the 
clusters initially fill a tidal radius of r t (R p ) following Baumgardt 
(1998), but also discuss the case where the clusters initially fill a 
tidal radius of r t (R a )- For the initial stellar mass function within 
each cluster, we adopt the model by Kroupa (2001) with a mass 
range of 0.08-15 Mq, which is realized by 15 discrete mass com- 
ponents in our FP model. 



4 BEST-FIT INITIAL MF AND RD 
4.1 Evolution of individual GCs 

First, we briefly discuss the effects of initial e and i on the evolu- 
tion of the individual clusters. Fig. [TJa) compares the clusters on 

3 It appears that such a continous update of rt with a realistic orbit calcu- 
lation is the first ever attempt for FP models. 
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values from both M and R histograms, which are constructed by 
using 1 1 bins between 10 3,75 M and 10 6,5 M for M and 11 bins 
between 10 2 ' 6 pc and 10 4 ' 8 pc for R, both equally spaced in loga- 
rithmic scale. Since the model M and R histograms are both con- 
strained by the observed number of clusters in our binning ranges, 
the degree of freedom for our % 2 test is 200 

Unlike M, R values evolve oscillating between R p and R a , 
and thus the model RD at 13 Gyr constructed from our 720 FP cal- 
culations may suffer a significant random noise. To decrease this 
noise, we build the model RD by adding up the probability distri- 
bution between R p and R a that are given by the orbital information 
at 13 Gyr. 



Figure 1. Evolution of cluster mass M relative to the initial mass Mi for 
a few clusters to show its dependence on the orbital eccentricity (a) and 
inclination (b). The values given in the plot are initial cluster parameters. 



eccentric orbits to those on circular orbits that have a radius of R a 
of the eccentric orbits, for two different M's. In case the initial r t 
is determined at R a , both massive and light clusters on eccentric 
orbits evolve faster than the cluster on a circular orbit, because the 
clusters on eccentric orbits are exposed to stronger tidal fields near 
R p . In case the initial r t is determined at R a , both massive and 
light clusters on eccentric orbits also evolve faster than the eccen- 
tric case where the initial r t is determined at R p , because for a 
given AI, the cluster with an initially larger r t would be more vul- 
nerable to the strong tidal field near R p . However, in case the initial 
rt is determined at R p , the light cluster evolves faster than the cir- 
cular case while the massive cluster evolves slower than the circular 
case. This is because the light cluster has a shorter relaxation time 
and thus its stars fill r t more quickly than the massive cluster when 
passing near R a . When the cluster passes near R p , the light cluster 
then loses stars more quickly over its shrunk r t . 

Fig. |TJb) shows that clusters with i > 45° evolve on similar 
timescales, but a cluster with a small i evolves somewhat faster due 
to a longer time spent while crossing the disc, which results in a 
stronger disc shock. 



4.2 Synthesis of FP calculations 

As discussed in ^3] we perform a total of 720 FP calcula- 
tions with different initial cluster conditions in 4-dimensional 
parameter space, M, R, e, and i. The goal of our present 
study is to find the best-fitting initial distributions of these vari- 
ables, and for this, we adopt a lognormal function dN(M) oc 
exp{— 0.5[log(M/Mp)/cri og M] 2 }dAf and a truncated power-law 
function diV(M) oc M~ a dM (only if M > Mi) for M and a 
softened power-law function dN(R) oc 47ri? 2 d7?/[l + (R/Rof] 
for R. We assume that the initial MF is independent of R. For the 
sake of simplicity, we do not parameterize the distributions for e 
and i, and adopt fixed isotropic distributions, diV(e) = 2ede and 
dN(i) = sin i di, respectively. 

Once the calculations of 720 FP models are done, the afore- 
mentioned set of initial MF and RD models are used to search for 
the best-fitting initial GC distributions in 4-dimensional parameter 
space: (M p , a\ 0& M, P, Ro) for the lognormal MF and (a, Mi, f3, 
Ro) for the power-law MF. We synthesize our 720 FP calculations 
with appropriate weights to produce a given initial MF and RD, and 
find a set of (M p , o"i og m, P, Ro) or (a, Mi, /3, Ro) that best fits the 
present-day MF and RD of the Galactic GC system. When finding 
the best-fitting initial MF and RD, we minimize the sum of the x 2 



4.3 Best-fitting initial GC distributions 

The best-fitting initial GC MFs and RDs from the \ 2 test between 
our models and the observed Galactic native (OH+BD) GCs are 
presented in TableQ] The best-fitting models for our standard initial 
condition have acceptably high p values (significances), and show 
a good agreement at 13 Gyr with the observed MF and RD (see 
Fig.|2j. Fig.|3ja) shows that the confidence intervals for 1, 2 and 3a 
are formed in a relatively small region, implying that our test gives 
rather small uncertainties. 

To see if our best-fitting models agree with the observed R de- 
pendence of the MF, we calculate \ 2 for the differences of (log M) 
and criog m between the model and the observation: 
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where subscripts o and m stand for the observation and the model, 
respectively, and subscript j represents the radial bin (we use thee 
radial bins for this test). We find that the p values from the above 
X 2 tests for our best-fitting lognormal and power-law models are 
larger than 50 per cent, implying that our best-fitting models have 
13 Gyr MFs whose R dependence is not significantly different from 
the observed MF (see Fig.[3p and c). 

Our best-fitting models have initial total masses in GCs (Mr,i) 
of 1.5-1.8 x 10 s M and the masses that have left GCs during the 
lifetime of the Galaxy (AM T ) of 1.2-1.5 x 10 8 M . Our AM T 
values are similar to or a few times larger than previous estimates 
by Baumgardt (1998; 4-9.5 x 10 7 M ) and Vesperini (1998; 5.5 x 
10 7 M ). Our best-fitting models result in the fraction of stars that 
remain in GCs until today of 14-18 per cent, which is comparable 
to the value obtained by Fall & Zhang (2001) for their lognormal 
model, 16 per cent. 

In spite of the comparable Mx,i, the MF of our best-fitting 
lognormal model is located in the more massive side with a nar- 
rower width than previous studies: it has a larger M p (10 M©) 



4 We have tried a 2D Kolmogorov-Smirnov test and a x 2 test with 2D 
bins as these tests can consider any correlations between M and R, but 
found that these 2D goodness-of-fit tests with a relatively small number of 
observed incidences, 95, result in a rather large acceptable ranges of param- 
eter space. Thus we first find the best-fitting parameters without considering 
the correlation between M and R instead, and then check if the best-fitting 
13 Gyr model MFs have the R dependence consistent with that of the ob- 
served MF. This way, we were able to find the model that reproduces the 
observed MF and RD very well simultaneously. 
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Table 1. Best-fitting parameters for initial GC distributions 



Model 

Standard 
Standard 
Circular 

cos i di 
r T {Ra) 



MF 

L 
P 
L 
L 
L 



log M„ 



°"log M 



log Mi 



13 



5.6r 

5.41; 
5.73; 

6.08" 



.11 

.15 



.11 
.18 
.10 
.0!) 
.06 
.00 



0.33 

0.48; 
0.29 
0.23 



-.00 
-.05 



.00 
.04 
+ .04 
.04 
+ .05 
.05 



2.31 



-.09 
-.09 



5.59 



+ .15 
18 



4.24; 

4.00; 
4.17; 

3.84; 
4.50" 



.21 
.21 
.22 
. 18 
.17 
.18 
.21 
.21 
.Hi 
.16 



R, 

2.9; 
2.0; 
2.0; 
2.0; 

2.7" 







M T ,i 

1.5 x 10 8 

1.8 x 10 8 

1.6 x 10 8 
1.6 x 10 8 

6.9 x 10 8 



X 

12.9 
14.8 
27.4 
16.2 
19.9 



p- value 

88 % 
79 % 
12 % 
70% 
46% 



Best-fitting initial distribution parameters for standard and non-standard initial conditions. 'L' in the MF column stands for the log- 
normal initial MF model, and 'P' for the power-law initial MF model. The p value is the probability of having a x 2 that is larger than 
the value obtained from our x 2 test between the model and the observation, whose degree of freedom is 20. The masses are in units 
of Mq and the radii in units of kpc. 
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Figure 2. M (upper panels) and R (lower panels) histograms of our best- 
fitting models at 13 Gyr (thick solid lines) with a lognormal initial MF (left- 
hand panels) and with a power-law initial MF (right-hand panels) along 
with the histograms for observed native GCs (thin solid lines). Also shown 
are the initial M and initial phased-mixed R histograms of our best-fitting 
models (dotted lines). 



and a smaller a losM (0.33) than Vesperini (1998; M p = 10 5 Mq, 
criogM =0.7) and Fall & Zhang (2001; M v = 2 x 1O 5 M , 
ciogM = 0.5). More various and realistic disruption mechanisms 
considered in our FP models increased the M of individual GCs 
and resulted in a larger initial M p to fit the present-day GCMF. 
Since our best-fitting model needs to fit the present-day GCMF 
with an increased M p , a smaller o"i og m is necessary as a compen- 
sation. 

Our best-fitting lognormal model has an initial RD of (3 — 
4.24, Rq = 2.9 kpc, but recall that this is the distribution of initial 
R a , and the true RD can be obtained by mixing the orbital phases of 
individual GCs. We find that the initial phase-mixed RD (see Fig.[2} 
is well described by a softened power-law function with (3 — 3.99, 
Ro — 1.7 kpc. The 13 Gyr RD of our best-fitting lognormal model 
in the figure shows a significant evolution from the initial value 
only in the small to intermediate R regions. This shows that the 
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log (M p /M s ) 
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Figure 3. Contours of A\ 2 values at 1, 4, and 9 (corresponding to 1, 2, 
and 3cr confidence intervals; solid lines) from the x 2 test between our log- 
normal model and the observation in the log Af p -<-ji og m plane (panel a). 
Also plotted are the 5 per cent significance levels for x 2 ({l°g A*0) (\ on S 
dash) and X^^logAf) (short dash), where the perpendicular tick marks 
point in the downhill direction. The minimum x 2 point (or the best-fitting 
model), marked with an asterisk, is located inside the 5 per cent significance 
levels of these two quantities, implying that our best-fitting model is con- 
sistent with the R, dependence of the observed MF. (log M) (panel b) and 
CTlog m values (panel c) of the best-fitting lognormal model (asterisks) and 
the observation (open circles) for three radial bins. Errorbars indicate the 
lcr uncertainties for the model. 



GCs with smaller R values are more vulnerable to mass loss and 
disruption because they have smaller relaxation times for a given 
M and encounter disc/bulge shocks more often. 

The MF of our best-fitting power-law model has a — 2.31, 
Mi = 10 5,59 Mq . This a value is in a good agreement with the ob- 
served range of a for giant molecular clouds and their star-forming 
cores in the local group of galaxies (e.g., Rosolowsky 2005), 1.5- 
2.5. Our best-fitting Mi is larger than that found by Parmentier 
& Gilmore (2005) for their truncated power-law IGCMF model, 
~ 10 Mq. We attribute this difference also to the fact that our 
models consider more various and realistic disruption mechanisms. 
When Mi < 10 4 ' 5 Mq, our 13 Gyr MF has either too small M v 
or too large ai og M, compared to the observed MF. Our best-fitting 
power-law model is found to have an initial phase-mixed RD with 
(3 — 3.91, Ro = 1.3 kpc, which is not too different from those of 
the lognormal model. 



4.4 Model dependences 

Now we discuss how the assumptions adopted in our models affect 
our results. In this subsection, we concentrate on our lognormal 
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Figure 4. GCMFs at 13 Gyr for several different initial conditions (lines) 
along with the observed GCMF (filled circles). Solid lines are our best- 
fitting lognormal GCMFs with our standard initial conditions. Each of the 
dashed and dotted lines shows the GCMF obtained with one different initial 
condition from the standard model as labeled. 



IGCMF models only. Fig.BJi shows that the 13 Gyr GCMF does not 
change significantly even when we have somewhat different initial 
distributions for e and i (circular orbits only instead of 2e de and 
cosidi instead of sinidz). The weak dependence of the GCMF 
on i distribution is easily expected from the evolution of individual 
clusters (Fig.Q}>), but the relatively weak dependence of the GCMF 
on e distribution is a bit surprising, considering the strong depen- 
dence on e seen in the M evolution of individual clusters (Fig.[T}i). 
The latter is thought to be because massive and light clusters have 
an opposite dependence of M on e (decreasing M with increasing 
e for massive clusters, but increasing M with increasing e for light 
clusters) and these opposite effects nearly cancel out each other. 

On the other hand, the evolution of the GCMF is quite sen- 
sitive to the choice of initial rt and Wo (see Fig.[4f>). As seen in 
Fig. [TJa), clusters evolve considerably faster when the inital r t is 
determined at R a than at R p , thus the 13 Gyr GCMF becomes 
significantly smaller in height in case of the former. The GCMF 
significantly shifts toward the lower M side also when the initial 
Wo is set to four instead of seven (when the cluster is initially less 
concentrated). This is because a cluster with a smaller Wo has rel- 
atively more stars in the outskirt compared to the core, and these 
stars have more chances of escaping the cluster than those in the 
core. 

Table [T] shows the best-fitting model parameters for our non- 
standard initial conditions as well. Mt,% is not sensitive to the ini- 
tial e or i distribution, but is more than four times larger when the 
initial rt is defined at R a . The mass that has left from GCs for 
the latter case, 6.6 x 10 s M , is a significant fraction of the cur- 
rent mass in the stellar halo between 4 and 25 kpc, 9 x 10 8 Mq 
(Suntzeff, Kinman, & Kraft 1990). This implies that the initial size 
of the GCs, along with the remnant gas expulsion in the early stage 
of the GC formation (Parmentier & Gilmore 2007), is an important 
factor in determining the origin of the halo stars. 

We were not able to find an initial MF and RD model that 
matches the observed Galactic GC system with a p value larger 
than 1 per cent, when the initial Wo is 4, instead of 7. This indicates 
that, within the limit of our parameterization for the initial MF and 
RD, a significant portion of the GCs are likely to have formed with 
Wo > 7. 



5 SUMMARY 

We have calculated the MF and RD evolution of the Galactic GC 
system using the most advanced and realistic FP model. By simul- 
taneously fitting both MF and RD of the observed GC system, we 
found that our best-fitting models have a higher M p for a lognormal 
initial MF and a higher Mi for a power-law initial MF than previous 
estimates, but our Mt/s are comparable to the previous results. 
Our best-fitting power-law MF model has a a = 2.2, which is in a 
good agreement with the observed range of a for giant molecular 
clouds and their star-forming cores. Our best-fitting lognormal and 
power-law models have initial phase-mixed RDs with ~ 4. Our 
best-fitting models, which are based on the assumptions of isotropic 
e and i distributions, agree well not only with the observed MF and 
RD, but also with the observed R dependence of the MF. Our re- 
sults are insensitive to the initial distributions of e and i, but are 
rather sensitive to how the initial tidal radius is defined and to the 
initial concentration of the clusters. If the clusters are assumed to 
be formed at the apocentre while filling the tidal radius there, ~ 75 
per cent of the mass in the current stellar halo could be attributed 
to the GCs as their origin. This implies that almost all of the mass 
in the stellar halo could be attributed to the GCs on eccentric orbits 
whose initial r t is much larger than r t (R p ), as well as to the GCs 
disrupted in the early stage of the cluster formation by the rem- 
nant gas expulsion. Finally, it appears that the clusters need to have 
a moderate-to-high initial concentration (Wo si 7) to explain the 
present-day MF and RD. 
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